#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stdbool.h>

#include "util.h"
//#include "ibm2d.h"
//#include "nsmg.h"
//#include "boundary.h"
//#include "boundary_list.h"
//#include "dlist.h"
//#include "output.h"

void multgf(double **f[],double **ar[],double **su[],double **res,int rf1,double **ae[],double **aw[],double **an[],double **as[],
            double **ap[], int *ni, int *nj, int *it, int *jt,double *residf,double *errf,double *work);
void coeff(int ih, double **ae, double **aw, double **an, double **as, double **ap, double **su, int ni, int nj, int it, int jt, double eps);
void restr1(double **phi1, int it1, int jt1, double **phi2, int ni2, int nj2, int it2, int jt2);
void restr2(double **phi1, int it1, int jt1, double **phi2, int ni2, int nj2, int it2, int jt2);
void interp(double **phi1, int it1, int jt1, double **phi2, int ni2, int nj2, int it2, int jt2);
void pgs(double **phi, double **su, double **ar, double rf1, int errf_i, double **ae, double **aw, double **an, double **as, double **ap,
           int ni, int nj, int it, int jt, int it1, int jt1, int residf_i,double *residf,double *errf);
void lgs(double **phi, double **su, double **ar, double rf1, int errf_i, double **ae, double **aw, double **an, double **as, double **ap,
           int ni, int nj, int it, int jt, int it1, int jt1, int residf_i,double *residf,double *errf);
void zlgs(double **phi, double **zsu, double **zar, double rf1, int errf_i, double **zae, double **zaw, double **zan, double **zas, double **zap,
           int ni, int nj, int it, int jt, int it1, int jt1, int residf_i,double *residf,double *errf);

//---subroutine
//multgf()  用于求解 Laplece 方程的多重网格方法 (FAS) 的驱动程序
//coeff()   该程序用于生成离散拉普拉斯方程的系数
//restr1()  此例程用于使用从较细网格到较粗网格的注入来重新 调整残差
//restr2()  此例程用于使用从较细级别网格到较粗级别网格的全权重重新 调整变量
//interp()  该例程用于通过使用从粗网格到细网格的双线性插值来执行 粗网格校正传输
//pgs()     该程序用于做点gauss-seidel松弛
//lgs()     该程序用于做线gauss-seidel松弛
//zlgs()    此例程用于执行斑马线gauss-seidel松弛
//print()

//FIRST EXECUTABLE STATEMENT LAP2D
int imax,jmax,km,itfas;
double tol,workm,rf1,eps;
int n[7];

int main()
{
    FILE *fp;
    if((fp=fopen("output1","w+")) == NULL)
    {
        printf("The file can not be opened.\n");
        exit(1);
    }
    //fprintf(fp,"it's work.");

    int it[7],jt[7];

    double **f[7];
    double **ae[7],**aw[7],**an[7],**as[7],**ap[7],**su[7],**ar[7];
    double **res;
    double *x,*y;
    
    int ni[7],nj[7];
    bool multig;
    char *hedfi[4] = {"fi ","pote","nt. ","fun."},
        *hedres[4] = {"res ","resi","dual","   "};

    double residf[7], errf[7], work=0.0;  //multgf(residf,errf,work)

    it[1]=65;
    jt[1]=65;
    it[2]=33;
    jt[2]=33;
    it[3]=17;
    jt[3]=17;
    it[4]=9;
    jt[4]=9;
    it[5]=5;
    jt[5]=5;
    it[6]=3;
    jt[6]=3;
//----------------------分配内存--------------------------------
    //double **f[7];
    for(int i=1; i<7;i++)
        f[i] = dmatrix(1, it[i], 1, jt[i]);  //注意释放内存；

    //double **ae[7],**aw[7],**an[7],**as[7],**ap[7],**su[7],**ar[7];
    //double **res;
    ae[1] = dmatrix(1, it[1], 1, jt[1]);  //注意释放内存；
    aw[1] = dmatrix(1, it[1], 1, jt[1]);
    an[1] = dmatrix(1, it[1], 1, jt[1]);
    as[1] = dmatrix(1, it[1], 1, jt[1]);
    ap[1] = dmatrix(1, it[1], 1, jt[1]);
    su[1] = dmatrix(1, it[1], 1, jt[1]);
    ar[1] = dmatrix(1, it[1], 1, jt[1]);
    res = dmatrix(1, it[1], 1, jt[1]);

    for(int i=2; i<7;i++)
    {
        ae[i] = dmatrix(1, it[i], 1, jt[i]);  //注意释放内存；
        aw[i] = dmatrix(1, it[i], 1, jt[i]);
        an[i] = dmatrix(1, it[i], 1, jt[i]);
        as[i] = dmatrix(1, it[i], 1, jt[i]);
        ap[i] = dmatrix(1, it[i], 1, jt[i]);
        su[i] = dmatrix(1, it[i], 1, jt[i]);
        ar[i] = dmatrix(1, it[i], 1, jt[i]);
    }

    //double *x,*y;
    x = dvector(1, it[1]);
    y = dvector(1, jt[1]);
//------------------------分配内存-------------------------------------

//--- default parameters
    multig = true;
    imax = 65;
    jmax = 65;
    km = 1;
    itfas = 10;
    workm = 2.0;
    if(multig == false)  workm = 0.0;
    eps = 1.0;
    rf1 = 1.0; // where to use ?

//--- control parameters display


//--- grid numbers for all levels
    ni[1] = imax;
    for(int i=2; i<7;i++)
    {
        ni[i] = (ni[i-1] - 1)/2 + 1;
    }

    nj[1] = jmax;
    for(int i=2; i<7;i++)
    {
        nj[i] = (nj[i-1] - 1)/2 + 1;
    }

    tol = 1e-8;

//--- number of relaxations at each level
    //n[7] = {1,1,1,1,1,1,1};
	int relax_n = 1;
	for (int i = 1; i <= 6; i++)
	{		
		n[i] = relax_n;
	}

//------coordinates  系数
    x[1]=1.0;
    for(int i=2; i<=imax; i++)
        x[i] = x[i-1] + 1.0;

    y[1]=1.0;
    for(int i=2; i<=imax; i++)
        y[i] = y[i-1] + 1.0;

//-----relaxation methods display  输出
    //if(km == 1)


//-----initial guess
    for(int i=1; i<=ni[1]; i++)
    {
        f[1][i][nj[1]] = 1.0;
        for(int j=1; j<=nj[1]-1; j++)
            f[1][i][j] = 0.0;
    }
    f[1][1][nj[1]] = 0;
    f[1][ni[1]][nj[1]] = 0;

//-----assembly coefficients (for linear problem)
    for(int ih = 1,i=1; ih<=32; ih=ih*2,i++)
    {
        coeff(ih,ae[i],aw[i],an[i],as[i],ap[i],su[i],ni[i],nj[i],it[i],jt[i],eps);
        //if( (ni[i+1]<3) || (nj[i+1]<3)) break;
    }

//-----Apply multigrid methods to solve finite difference equations
    multgf( f, ar, su, res, rf1, ae, aw, an, as,
             ap,  ni,  nj,  it,  jt, residf, errf, &work);
    printf(" imax=%d  jmax=%d   work=%lf  residf=%.10lf  errf=%.10lf  \n\n", 
                imax, jmax, work, residf[1], errf[1]);

    double resido, rate, workt=0.0;
    for(int ite=1; ite<=itfas; ite++)
    {
        resido = residf[1];
        multgf( f, ar, su, res, rf1, ae, aw, an, as,
               ap, ni, nj, it,  jt, residf, errf, &work);
        rate = residf[1]/resido;
        workt = workt + work;
        printf("ite=%d  imax=%d  jmax=%d  workt=%lf  work=%lf  residf=%.10lf  errf=%.10lf  rate=%lf\n\n", 
                ite, imax, jmax, workt,work,residf[1],errf[1],rate);
        if(residf[1] < tol) break;
    }
    //printf("imax=%d  jmax=%d  workt=%lf  work=%lf  residf=%lf  errf=%lf  rate=%lf\n", 
                //imax, jmax, workt,work,residf[1],errf[1],rate);

//-----final output - potential function and the residuals  输出

//format statement   输出格式声明

    fclose(fp);

//----------------------释放内存--------------------------------
    //double **f[7];
    for(int i=1; i<7;i++)
    free_dmatrix(f[i], 1, it[i], 1, jt[i]);

    //double **ae[7],**aw[7],**an[7],**as[7],**ap[7],**su[7],**ar[7];
    //double **res;
    free_dmatrix(ae[1], 1, it[1], 1, jt[1]);
    free_dmatrix(aw[1], 1, it[1], 1, jt[1]);
    free_dmatrix(an[1], 1, it[1], 1, jt[1]);
    free_dmatrix(as[1], 1, it[1], 1, jt[1]);
    free_dmatrix(ap[1], 1, it[1], 1, jt[1]);
    free_dmatrix(su[1], 1, it[1], 1, jt[1]);
    free_dmatrix(ar[1], 1, it[1], 1, jt[1]);
    free_dmatrix(res, 1, it[1], 1, jt[1]);

    for(int i=2; i<7;i++)
    {
        free_dmatrix(ae[i], 1, it[i], 1, jt[i]);
        free_dmatrix(aw[i], 1, it[i], 1, jt[i]);
        free_dmatrix(an[i], 1, it[i], 1, jt[i]);
        free_dmatrix(as[i], 1, it[i], 1, jt[i]);
        free_dmatrix(ap[i], 1, it[i], 1, jt[i]);
        free_dmatrix(su[i], 1, it[i], 1, jt[i]);
        free_dmatrix(ar[i], 1, it[i], 1, jt[i]);
    }

    //double *x,*y;
    free_dvector(x, 1, it[1]);
    free_dvector(y, 1, jt[1]);
//------------------------释放内存-------------------------------------
    
    printf("finished!!\n");

}

//double multgf(double residf[], double errf, double work)
void multgf(double **f[],double **ar[],double **su[],double **res,int rf1,double **ae[],double **aw[],double **an[],double **as[],
            double **ap[], int *ni, int *nj, int *it, int *jt,double *residf,double *errf,double *work)
{
    //***FIRST EXECUTABLE STATEMENT MULTGF
    printf("1111");
    *work = 0.0;
    int flag = 0;
    do{
        //-----relaxation on grid level 1
        int level = 1;
        int ih=1;
        //coeff(ih,ae[level],aw[level],an[level],as[level],ap[level],su[level],ni[level],nj[level],it[level],jt[level],eps);
        for(int i=1; i<=n[1]; i++) 
        {
            if(km == 1) pgs(f[1],su[1],res,rf1,1,ae[1],aw[1],an[1],as[1],ap[1],ni[1],nj[1],it[1],jt[1],it[1],jt[1],1,residf,errf);
            if(km == 2) lgs(f[1],su[1],res,rf1,1,ae[1],aw[1],an[1],as[1],ap[1],ni[1],nj[1],it[1],jt[1],it[1],jt[1],1,residf,errf);
            if(km == 3) zlgs(f[1],su[1],res,rf1,1,ae[1],aw[1],an[1],as[1],ap[1],ni[1],nj[1],it[1],jt[1],it[1],jt[1],1,residf,errf);
            *work = *work + 1.0;
            printf("residph=%.10lf \n",residf[1]);
            //printf("work1=%lf\n",work);
        }

        //printf("level=%d work=%lf  residf=%lf \n", level,work,residf[1]);

        //grid level 2-6
        for(level=2; level<=6; level++)
        {
            ih = pow(2,level-1);
            //printf("level=%d work=%lf  residf=%lf \n", level,work,residf[1]);
            //printf("coeff\n");
            //coeff(ih,ae[level],aw[level],an[level],as[level],ap[level],su[level],ni[level],nj[level],it[level],jt[level],eps);

            //-----restriction
            restr1(res,it[1],jt[1],ar[level],ni[level],nj[level],it[level],jt[level]);
            restr2(f[level-1],it[level-1],jt[level-1],f[level],ni[level],nj[level],it[level],jt[level]);

            //--- coarse grid right hand side
            for(int i=2; i<=ni[level]-1; i++)
                for(int j=2; j<=nj[level]-1; j++)
                {
                    ar[level][i][j] = ar[level][i][j] - (ae[level][i][j]*f[level][i+1][j] + aw[level][i][j]*f[level][i-1][j]
                                    + an[level][i][j]*f[level][i][j+1] + as[level][i][j]*f[level][i][j-1] - ap[level][i][j]*f[level][i][j]);
                }

            //-----relaxation on grid level
            for(int i=1; i<=n[level]; i++)
            {
                if(km == 1) pgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                if(km == 2) lgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                if(km == 3) zlgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                *work = *work + 1.0/pow(4.0,level-1);
                printf("residph=%.10lf \n",residf[level]);

            }

            if(ni[level]<3 || nj[level]<3)  {  flag=1; break;}
        }
        //printf("2-6work=%lf\n",work);
        //printf("level = %d ",level);
        
/*
//--------------------------------------------------------------------------------------
        printf("level = %d ",level);
        if(flag == 0)
        {
        //level++;
        ih = pow(2,level-1);
            //printf("level=%d work=%lf  residf=%lf \n", level,work,residf[1]);

            //-----restriction
            restr1(res,it[1],jt[1],ar[level],ni[level],nj[level],it[level],jt[level]);
            restr2(f[level-1],it[level-1],jt[level-1],f[level],ni[level],nj[level],it[level],jt[level]);

            //--- coarse grid right hand side
            for(int i=2; i<=ni[level]-1; i++)
                for(int j=2; j<=nj[level]-1; j++)
                {
                    ar[level][i][j] = ar[level][i][j] - (ae[level][i][j]*f[level][i+1][j] + aw[level][i][j]*f[level][i-1][j]
                                    + an[level][i][j]*f[level][i][j+1] + as[level][i][j]*f[level][i][j-1] - ap[level][i][j]*f[level][i][j]);
                }

            //-----relaxation on grid level
            for(int i=1; i<=n[level]; i++)
            {
                if(km == 1) pgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level);
                if(km == 2) lgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level);
                if(km == 3) zlgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level);
                work = work + 1.0/pow(4.0,level-1);
                //printf("residph=%.10lf \n",residf[level]);
            }

            if(ni[level]<5 || nj[level]<3)  {  flag=1; break;}
        }
        //----------------------------------------------------------------
        if(flag == 0)
        {
            level++;
            printf("level = %d ",level);
            ih = pow(2,level-1);
            //printf("level=%d work=%lf  residf=%lf \n", level,work,residf[1]);

            //-----restriction
            restr1(res,it[1],jt[1],ar[level],ni[level],nj[level],it[level],jt[level]);
            restr2(f[level-1],it[level-1],jt[level-1],f[level],ni[level],nj[level],it[level],jt[level]);

            //--- coarse grid right hand side
            for(int i=2; i<=ni[level]-1; i++)
                for(int j=2; j<=nj[level]-1; j++)
                {
                    ar[level][i][j] = ar[level][i][j] - (ae[level][i][j]*f[level][i+1][j] + aw[level][i][j]*f[level][i-1][j]
                                    + an[level][i][j]*f[level][i][j+1] + as[level][i][j]*f[level][i][j-1] - ap[level][i][j]*f[level][i][j]);
                }

            //-----relaxation on grid level
            for(int i=1; i<=n[level]; i++)
            {
                if(km == 1) pgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level);
                if(km == 2) lgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level);
                if(km == 3) zlgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level);
                work = work + 1.0/pow(4.0,level-1);
                //printf("residph=%.10lf \n",residf[level]);
            }

            if(ni[level]<3 || nj[level]<3)  break;
        }
        level++;*/
//---------------------------------------------------------------------------------------
        //printf("level = %d   ",level);
        level--;
        //printf("level = %d\n",level);
        if( level == 6) //level 6
        {
            coeff(ih,ae[level],aw[level],an[level],as[level],ap[level],su[level],ni[level],nj[level],it[level],jt[level],eps);
            restr2(f[level-1],it[level-1],jt[level-1],su[level],ni[level],nj[level],it[level],jt[level]);
            for(int i=1; i<ni[level]; i++)
                for(int j=1; j<nj[level]; j++)
                {
                    su[level][i][j] = f[level][i][j] - su[level][i][j];
                }
            interp(f[level-1],it[level-1],jt[level-1],su[level],ni[level],nj[level],it[level],jt[level]);
            level--;
        }

        for(; level>1; level--)  //level 2-5
        {
            //relaxation
            for(int i=1; i<=n[level]; i++)
            {
                if(km == 1) pgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                if(km == 2) lgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                if(km == 3) zlgs(f[level],ar[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                *work = *work + 1.0/pow(4.0,level-1);
                printf("residph=%.10lf \n",residf[level]);
            }

            //-----coarse grid correction
            coeff(ih,ae[level],aw[level],an[level],as[level],ap[level],su[level],ni[level],nj[level],it[level],jt[level],eps);
            restr2(f[level-1],it[level-1],jt[level-1],su[level],ni[level],nj[level],it[level],jt[level]);
            for(int i=1; i<ni[level]; i++)
                for(int j=1; j<nj[level]; j++)
                {
                    su[level][i][j] = f[level][i][j] - su[level][i][j];
                }
            interp(f[level-1],it[level-1],jt[level-1],su[level],ni[level],nj[level],it[level],jt[level]);
        }
        //printf("1---level = %d\n",level);
        if( level == 1) //level 1
        {
            ih=1;
            coeff(ih,ae[level],aw[level],an[level],as[level],ap[level],su[level],ni[level],nj[level],it[level],jt[level],eps);
            for(int i=1; i<=n[level]; i++)
            {
                if(km == 1) pgs(f[level],su[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                if(km == 2) lgs(f[level],su[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                if(km == 3) zlgs(f[level],su[level],res,rf1,level,ae[level],aw[level],an[level],as[level],ap[level],ni[level],nj[level],it[level],jt[level],it[1],jt[1],level,residf,errf);
                *work = *work + 1.0/pow(4.0,level-1);
                printf("residph=%.10lf \n",residf[1]);
            }
        }
        if(*work>workm) {printf("mult_work=%lf\n\n",*work); break;}
    }while(  (*work<=workm) || (fabs(residf[1])>tol));

    //return residf[1];
}

//coeff(ih, ae, aw, an, as, ap, su, ni, nj, it, jt, eps)
void coeff(int ih, double **cae, double **caw, double **can, double **cas, double **cap, double **csu, int ni, int nj, int it, int jt, double eps)
{
    float hi = (float)ih;
    float h2r = 1.0/hi/hi;
    for(int i=2; i<=ni-1; i++)
        for(int j=2; j<=nj-1; j++)
        {
            cae[i][j] = 1.0*eps*h2r;
            caw[i][j] = 1.0*eps*h2r;
            can[i][j] = 1.0*h2r;
            cas[i][j] = 1.0*h2r;
            cap[i][j] = cae[i][j] + caw[i][j] + can[i][j] + cas[i][j];
            csu[i][j] = 0.0;
        }
}

//restr1(phi1,it1,jt1,phi2,ni2,nj2,it2,jt2)
void restr1(double **phi1, int it1, int jt1, double **phi2, int ni2, int nj2, int it2, int jt2)
{
    for(int i=2; i<=ni2-1; i++)
    {
        int i1 = 2*(i-1) + 1;
        //int work=0;
        for(int j=2; j<=nj2-1; j++)
        {
            //work++;
            int j1 = 2*(j-1) - 1;
            phi2[i][j] = 0.25*phi1[i1][j1] + 0.125*(phi1[i1+1][j1] + phi1[i1-1][j1] - phi1[i1][j1+1] + phi1[i1][j1-1])
                       + 0.0625*(phi1[i1+1][j1+1] - phi1[i1-1][j1+1] + phi1[i1+1][j1-1] + phi1[i1-1][j1-1]);
        }
    }

}

//restr2(phi1,it1,jt1,phi2,ni2,nj2,it2,jt2)
void restr2(double **phi1, int it1, int jt1, double **phi2, int ni2, int nj2, int it2, int jt2)
{
    //calculate the coarser grid values for the boundary points
    for(int j=1; j<=nj2; j=j+nj2-1)
    {
        int j1 = 2*(j-1) + 1;
        for(int i=1; i<=ni2; i++)
        {
            int i1 = 2*(i-1) + 1;
            phi2[i][j] = phi1[i1][j1];
        }
    }

    for(int i=1; i<=ni2; i=i+ni2-1)
    {
        int i1 = 2*(i-1) + 1;
        for(int j=1; j<=nj2; j++)
        {
            int j1 = 2*(j-1) + 1;
            phi2[i][j] = phi1[i1][j1];
        }
    }

    //calculate the coarser grid values for the interior points
    for(int j=2; j<=nj2-1; j++)
    {
        int j1 = 2*(j-1) + 1;
        for(int i=2; i<=ni2-1; i++)
        {
            int i1 = 2*(i-1) + 1;
            phi2[i][j] = 0.5*phi1[i1][j1] + 0.125*(phi1[i1+1][j1] + phi1[i1-1][j1] + phi1[i1][j1+1] + phi1[i1][j1-1]);
        }
    }
}

//interp(phi1,it1,jt1,phi2,ni2,nj2,it2,jt2)
void interp(double **phi1, int it1, int jt1, double **phi2, int ni2, int nj2, int it2, int jt2)
{
    int ni1 = (ni2-1)*2 + 1;
    int nj1 = (nj2-1)*2 + 1;

    for(int i=3; i<=ni1-2; i=i+2)
    {
        int i2 = (i+1)/2;
        for(int j=3; j<=nj1-2; j=j+2)
        {
            int j2 = (j+1)/2;
            phi1[i][j] = phi1[i][j] + phi2[i2][j2];
        }
        for(int j=2; j<=nj1-1; j=j+2)
        {
            int j2 = (j+2)/2;
            phi1[i][j] = phi1[i][j] + 0.5*(phi2[i2][j2-1] + phi2[i2][j2]);
        }
    }

    for(int i=2; i<=ni1-1; i=i+2)
    {
        int i2 = (i+2)/2;
        for(int j=3; j<=nj1-2; j=j+2)
        {
            int j2 = (j+1)/2;
            phi1[i][j] = phi1[i][j] + 0.5*(phi2[i2-1][j2] + phi2[i2][j2]);
        }
        for(int j=2; j<=nj1-1; j=j+2)
        {
            int j2 = (j+2)/2;
            phi1[i][j] = phi1[i][j] + 0.25*(phi2[i2-1][j2-1] + phi2[i2-1][j2] + phi2[i2][j2-1] + phi2[i2][j2]);
        }
    }
}


//pgs(phi,su,ar,rf1,errph,ae,aw,as,ap,ni,nj,it,jt,it1,jt1,residph)  //errhp residph删掉
void pgs(double **phi, double **psu, double **par, double rf1, int errf_i, double **pae, double **paw, double **pan, double **pas, double **pap,
           int ni, int nj, int it, int jt, int it1, int jt1, int residf_i,double *residf,double *errf)
{
    residf[residf_i] = 0.0;
    errf[errf_i] = 0.0;
    //double residph = 0.0;
    //double errph = 0.0;

    //calculate the new phi(i,j)
    double phitem;
    double phinew;
    double phie;
    for(int i=2; i<=ni-1; i++)
        for(int j=2; j<=nj-1; j++)
        {
            phitem = psu[i][i] + pae[i][j]*phi[i+1][j] + paw[i][j]*phi[i-1][j] + pan[i][j]*phi[i][j+1] + pas[i][j]*phi[i][j-1];
            phitem = phitem/pap[i][j];
            phinew = rf1*phitem + (1.0-rf1)*phi[i][j];
            phie = phinew - phi[i][j];
            phi[i][j] = phinew;
            errf[errf_i] = (errf[errf_i]>fabs(phie)) ? (errf[errf_i]):(fabs(phie));
            //errph = dmax1(errph,dabs(phie));   ///dmax1///
        }

    //calculate the new residuals
    for(int i=2; i<=ni-1; i++)
        for(int j=2; j<=nj-1; j++)
        {
            par[i][j] = psu[i][j] + pae[i][j]*phi[i+1][j] + paw[i][j]*phi[i-1][j]
                        + pan[i][j]*phi[i][j+1] + pas[i][j]*phi[i][j-1] - pap[i][j]*phi[i][j];
            residf[residf_i] = residf[residf_i] + par[i][j]*par[i][j];
        }

    residf[residf_i] = sqrt(residf[residf_i]/(float)((ni-2)*(nj-2)));
    //residph = sqrt(residph/(float)((ni-2)*(nj-2)));
    //return residph;
}


//lgs(phi,su,ar,rf1,errph,ae,aw,as,ap,
//      ni,nj,it,jt,it1,jt1,residph)   //errhp residph删掉
void lgs(double **phi, double **lsu, double **lar, double rf1, int errf_i, double **lae, double **law, double **lan, double **las, double **lap,
         int ni, int nj, int it, int jt, int it1, int jt1, int residf_i,double *residf,double *errf)
{
    double a[66],b[66],c[66],d[66];

    residf[residf_i] = 0.0;
    errf[errf_i] = 0.0;

    //start w-e sweeps
    for(int i=2; i<=ni-1; i++)
    {
        c[1] = phi[i][1];

        //start s-n traverse
        int j=2;
        a[j] = lan[i][j];
        b[j] = las[i][j];
        c[j] = lae[i][j]*phi[i+1][j] + law[i][j]*phi[i-1][j] + lsu[i][j];
        c[j] = c[j] + b[j]*c[j-1];
        d[j] = lap[i][j];

        //calculate coefficients of recurrence formula
        double term = 1.0/d[j];
        a[j] = a[j]*term;
        c[j] = c[j]*term;

        for(int j=3; j<=nj-1; j++)
        {
            //assembly coefficients
            a[j] = lan[i][j];
            b[j] = las[i][j];
            c[j] = lae[i][j]*phi[i+1][j] + law[i][j]*phi[i-1][j] + lsu[i][j];
            d[j] = lap[i][j];

            //calculate coefficients of recurrence formula
            term = 1.0/(d[j] - b[j]*a[j-1]);
            a[j] = a[j]*term;
            c[j] = (c[j] + b[j]*c[j-1])*term;
        }

        //obtain new phi and errors
        double phitem;
        double phinew;
        double phie;
        for(int jj=2; jj<=nj-1; jj++)
        {
            int j = nj+1-jj;
            phitem = a[j]*phi[i][j+1] + c[j];
            phinew = rf1*phitem + (1.0-rf1)*phi[i][j];
            phie = phinew - phi[i][j];
            phi[i][j] = phinew;
            errf[errf_i] = (errf[errf_i]>fabs(phie)) ? (errf[errf_i]):(fabs(phie));
            //errph = dmax1(errph,dabs(phie));
        }
    }

    //calculate the new residuals
    for(int i=2; i<=ni-1; i++)
        for(int j=2; j<=nj-1; j++)
        {
            lar[i][j] = lsu[i][j] + lae[i][j] *phi[i+1][j] + law[i][j]*phi[i-1][j]
                       + lan[i][j]*phi[i][j+1] + las[i][j]*phi[i][j-1] - lap[i][j]*phi[i][j];
            residf[residf_i] = residf[residf_i] + lar[i][j]*lar[i][j];
        }

    residf[residf_i] = sqrt(residf[residf_i]/(float)((ni-2)*(nj-2)));
    //return residph;
}


//zlgs(phi,su,ar,rf1,errph,ae,aw,as,ap,
//      ni,nj,it,jt,it1,jt1,residph)  //errhp residph删掉
void zlgs(double **phi, double **zsu, double **zar, double rf1, int errf_i, double **zae, double **zaw, double **zan, double **zas, double **zap,
         int ni, int nj, int it, int jt, int it1, int jt1, int residf_i,double *residf,double *errf)
{
    double a[66],b[66],c[66],d[66];

    residf[residf_i] = 0.0;
    errf[errf_i] = 0.0;
    //double residph = 0.0;
    //double errph = 0.0;
    //printf("residph = %.10lf  errph = %.10lf\n", residf[residf_i] ,errf[errf_i]);

    double phitem;
    double phinew;
    double phie;
    double term;

    //start w-e along even number lines
    for(int i=2; i<=ni-1; i=i+2)
    {
        c[1] = phi[i][1];

        //start s-n traverse
        int j=2;
        a[j] = zan[i][j];
        b[j] = zas[i][j];
        c[j] = zae[i][j]*phi[i+1][j] + zaw[i][j]*phi[i-1][j] + zsu[i][j];
        c[j] = c[j] + b[j]*c[j-1];
        d[j] = zap[i][j];

        //calculate coefficients of recurrence formula
        term = 1.0/d[j];
        a[j] = a[j]*term;
        c[j] = c[j]*term;

        for( j=3; j<=nj-1; j++)
        {
            //assembly coefficients
            a[j] = zan[i][j];
            b[j] = zas[i][j];
            c[j] = zae[i][j]*phi[i+1][j] + zaw[i][j]*phi[i-1][j] + zsu[i][j];
            d[j] = zap[i][j];

            //calculate coefficients of recurrence formula
            term = 1.0/(d[j]-b[j]*a[j-1]);
            a[j] = a[j]*term;
            c[j] = (c[j] + b[j]*c[j-1])*term;
        }

        //obtain new phi and errors
        for(int jj=2; jj<=nj-1; jj++)
        {
            int j = nj+1-jj;
            phitem = a[j]*phi[i][j+1] + c[j];
            phinew = rf1*phitem + (1.0-rf1)*phi[i][j];
            phie = phinew - phi[i][j];
            phi[i][j] = phinew;
            errf[errf_i] = (errf[errf_i]>fabs(phie)) ? (errf[errf_i]):(fabs(phie));
            //errph = dmax1(errph,dabs(phie));
        }
    }

    //start w-e sweeps along odd number lines
    for(int i=3; i<=ni-2; i=i+2)
    {
        c[1] = phi[i][1];

        //start s-n traverse
        int j=2;
        a[j] = zan[i][j];
        b[j] = zas[i][j];
        c[j] = zae[i][j]*phi[i+1][j] + zaw[i][j]*phi[i-1][j] + zsu[i][j];
        c[j] = c[j] + b[j]*c[j-1];
        d[j] = zap[i][j];

        //calculate coefficients of recurrence formula
        term = 1.0/d[j];
        a[j] = a[j]*term;
        c[j] = c[j]*term;

        for( j=3; j<=nj-1; j++)
        {
            //assembly coefficients
            a[j] = zan[i][j];
            b[j] = zas[i][j];
            c[j] = zae[i][j]*phi[i+1][j] + zaw[i][j]*phi[i-1][j] + zsu[i][j];
            d[j] = zap[i][j];

            //calculate coefficients of recurrence formula
            term = 1.0/(d[j]-b[j]*a[j-1]);
            a[j] = a[j]*term;
            c[j] = (c[j] + b[j]*c[j-1])*term;
        }

        //obtain new phi and errors
        for(int jj=2; jj<=nj-1; jj++)
        {
            int j = nj+1-jj;
            phitem = a[j]*phi[i][j+1] + c[j];
            phinew = rf1*phitem + (1.0-rf1)*phi[i][j];
            phie = phinew - phi[i][j];
            phi[i][j] = phinew;
            errf[errf_i] = (errf[errf_i]>fabs(phie)) ? (errf[errf_i]):(fabs(phie));
            //errph = dmax1(errph,dabs(phie));
        }
    }


    //calculate the new residuals
    for(int i=2; i<=ni-1; i++)
        for(int j=2; j<=nj-1; j++)
        {
            zar[i][j] = zsu[i][j] + zae[i][j]*phi[i+1][j] + zaw[i][j]*phi[i-1][j]
                       + zan[i][j]*phi[i][j+1] + zas[i][j]*phi[i][j-1] - zap[i][j]*phi[i][j];
            residf[residf_i] = residf[residf_i] + zar[i][j]*zar[i][j];
        }

    residf[residf_i] = sqrt(residf[residf_i]/(float)((ni-2)*(nj-2)));

    //printf("residph = %.10lf  errph = %.10lf\n\n", residf[residf_i] ,errf[errf_i]);
    //printf("residph=%lf \n",residph);
    //return residph;
}

